c2.pps=function(data,ybari,M,param){
  data=as.name(data);
  ybari=ybari; n=length(ybari)
  if(param=='mean'){
    est=(1/n)*sum(ybari)
    vhat=(1/(n*(n-1)))*sum((ybari-est)^2)
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } else{
    est=(M/n)*sum(ybari)
    vhat=M^2*((1/(n*(n-1)))*sum((ybari-est)^2))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } 
  cat("","\n","Results from 2-stage Cluster sample (pps): Data =",data,"\n","\n","n =",n,"\n","Estimate of",param,"=",est,
  "\n","Vhat(",param,") =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",lower,"Upper Bound =",upper,"\n","")
  results=list(est=est,data=data,n=n,M=M,vhat=vhat,bound=bound,lower=lower,upper=upper,param=param)
}

# to use the function with its call:
# c2.pps(data,ybari,M,param)

# data: name of dataset, in quotes
# y: numeric vector; the total of all observations in the ith cluster
# m: numeric vector; number of elements in cluster i, i=1,2,...,N
# M: number of elements in population (single number)
# param: c('mean','total')